Binary Black Hole Mergers in 3d Numerical Relativity 
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The standard approach to the numerical evolution of black hole data using the ADM formulation 
with maximal slicing and vanishing shift is extended to non-symmetric black hole data containing 
black holes with linear momentum and spin by using a time-independent conformal rescaling based 
on the puncture representation of the black holes. We give an example for a concrete three dimen- 
sional numerical implementation. The main result of the simulations is that this approach allows 
for the first time to evolve through a brief period of the merger phase of the black hole inspiral. 
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I. INTRODUCTION 

While the binary black hole problem is essentially un- 
solved, there have been many advances in our under- 
standing of the general relativistic evolution problem. 
The situation of interest is how two black holes orbit each 
other, spiral inwards, eventually merge, and emit grav- 
itational waves in the process. Gravitational wave de- 
tectors are being built that eventually will detect waves 
from such sources M. Although at least as interesting 
from an astrophysical as opposed to purely general rel- 
ativistic view point, here we will not consider neutron 
stars or other matter sources (e.g. Q). To name just 
three areas of research related to the inspiral of two black 
holes: post-Newtonian methods are applicable if the sys- 
tem is not too strongly general relativistic ||, that is in 
the early phase of the inspiral; perturbation methods can 
model the final stages of the merger Q| ; and full numeri- 
cal relativity can be applied to the late strong-interaction 
phase if one imposes axisymmetry |||| . 

What may appear surprising is that the general, three- 
dimensional two black hole problem has not been solved 
numerically to some degree. After all, there has been 
a lot of theoretical work on the Einstein equations, and 
even though they are complicated, one might expect that 
eventually computers become powerful enough to at least 
allow some rough investigations into the general problem. 

However, general relativity has at least two charac- 
teristic features, gauge freedom and singularity forma- 
tion, which would prevent us from solving the problem 
completely with current computer codes even if an in- 
finitely fast computer became available. Here gauge free- 
dom refers to the freedom to choose coordinates, which 
in numerical relativity is a hard problem since in general 
one cannot specify coordinates for the entire spacetime 
in advance. For an adequate representation of the space- 
time on a numerical grid one usually has to construct the 
coordinates numerically as one of the steps in the numer- 



ical evolution scheme. In this dynamic construction of 
coordinates one has to avoid the formation of coordinate 
singularities, but even if these are absent, it is still possi- 
ble that completely regular initial data develops physical 
singularities, as is the case for typical black hole space- 
times. Both problems, the choice of coordinates and the 
representation of black hole spacetimes on a numerical 
grid, are still awaiting a completely satisfactory solution. 

What has been achieved in 3d numerical relativity is 
the evolution of weak gravitational waves (non- linear 
waves that are weak enough to not collapse to black 
holes), and the evolution of a single spherically sym- 
metric black hole ||-[nj in Cartesian coordinates and for 
non-trivial slicings. Recently, long term stable evolutions 
have been achieved for single black holes u^L A sin- 
gle moving black hole has been simulated in ||l3|, which 
should become important for black hole collisions. Ax- 
isymmetric collisions of two black holes with a 3d code 
have been reported in pj ] . Non-axisymmetric spacetimes 
containing a single distorted black hole are studied in 
|]l5f . Various other 3d projects are actively pursued, most 
notably by the US binary black hole alliance Jl6[ . 

In this paper we introduce a new approach that makes 
the study of non-axisymmetric binary black hole space- 
times possible, at least for short time intervals. The cru- 
cial new insight is that the general binary black hole ini- 
tial data of |17| can be evolved on R 3 without special 
inner boundary using a conformal rescaling that is con- 
stant in time. The main result is the evolution of the 
apparent horizon, which for appropriate black hole ini- 
tial data has two components initially that are replaced 
by a single component during evolution. 

The main limitation of the current implementation is 
that only a brief interval of the merger phase of the in- 
spiral of two black holes can be studied. This is due 
to the well-known problem of grid-stretching associated 
with our particular gauge choice, maximal slicing and 
vanishing shift. We outline below how the gauge condi- 
tions can be generalized. Here we want to demonstrate 
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our approach with one minimal example for general black 
hole data, a more detailed study is left to a future pub- 
lication. 

To summarize, initial data for two black holes of 
unequal mass that are moving and spinning are con- 
structed following |17[| . The evolution is carried out in 
the 3+1 ADM formulation with maximal slicing and van- 
ishing shift using BAM, a bifunctional adaptive mesh 
code for coupled elliptic and hyperbolic systems (only 
a fixed mesh refinement of nested grids is used here; see 
|H|[f7||l8) about earlier versions of this code). An evolu- 
tion time of about 7-10M (in units of the ADM mass) 
is achieved. We can study the geometric information in 
the collapse of the lapse function, and we find the appar- 
ent horizon and show data for the transition from two 
outermost marginally trapped closed surfaces to a single 
one. 

In the following, we describe the method and its imple- 
mentation, present the physical data, and conclude with 
a discussion. 



II. DESCRIPTION OF THE PUNCTURE 
METHOD 

We begin with the (3+l)-dimensional Arnowitt-Deser- 
Misner (ADM) decomposition of the 4-dimensional Ein- 
stein equations |fL9[ . The dynamical fields are the 3- 
metric g ab and its extrinsic curvature K ab on a 3-manifold 
E, both depending on space (points in S) and a time pa- 
rameter, t. The foliation of the 4-dimcnsional spacetime 
into hypersurfaces £ is characterized in the usual way by 
a lapse function a and a shift vector (3 a . The Einstein 
equations become 

d t 9ab = -2aK ab + D a /3 b + D b f3 a , (1) 
d t K ab = -D a D b a + a(R ab - 2K ac K c b + K ab K) 

+P c D c K ab + K ac D b /3 c + K cb D a f3 c , (2) 

= D b (K ab - g ab K), (3) 

= R-K ab K ab + K 2 , (4) 

where R ab is the 3-Ricci tensor, R the Ricci scalar, K 
the trace of the extrinsic curvature, and D a the covariant 
derivative derived from the 3-metric. 

We solve the constraints on an initial hypersurface at 
t = for conformally flat two black hole data by the 
method introduced in p"7| . The solution is found in terms 
of conformally transformed quantities, g ab = ip 4 5 ab and 
K ab = ip~ 2 K ab . On a Cartesian grid, one chooses two 
points el and c<i that represent the internal asymptoti- 
cally flat regions of the two black holes. 

The diffeomorphism constraint is solved by the Bowen- 
York extrinsic curvature |po[] with parameters P\ and Pi 
for the linear momenta and S\ and Si for the spins, 

2 3 

Kps = ^([^(P"^ + P b n a - (S ab - n a n b )P c n c )] z 



+ [^(e acd S c n d n b + e bcd S c n d n a )]i), (5) 

where n a is the radial normal vector in Cartesian coor- 
dinates, n a = x a /r, and where [■ ■ -]i denotes the term 
for the i-th black hole with i\ — \x — <3|. The van- 
ishing of the Hamiltonian constraint H, (^), becomes a 
non-linear elliptic equation for the conformal factor 
The novel feature of Jl7j is that the equation for ip is 
rewritten for i? 3 , where the points C\ and c% are part of 
the computational domain. The solution has the form 
tp = u + J2i=i m i/(2ri), where u is an everywhere regu- 
lar function. Recall that Schwarzschild initial data can 
be written as K ab = and ip = 1 4- m/(2r). 

We presented the initial data in some detail because 
this particular construction leads to a natural treatment 
of the inner boundary during evolution. On the initial 
slice the metric diverges at c[ and 62 (in the right man- 
ner to represent the inner asymptotically flat regions). 
As already done in f|,[l(|, an evolution of Schwarzschild 
initial data can be performed on i? 3 , if the metric is di- 
vided by the time-independent conformal factor ip A = 
(1 + m/(2r)) 4 , g ab = ip 4 g ab - The derivatives of the phys- 
ical metric g ab are obtained as the numerical derivatives 
of the regular conformal metric g ab plus terms contain- 
ing the derivatives of the conformal factor ip, which are 
computed analytically. Furthermore, the grid is offset so 
that the puncture is not one of the grid points. 

For binary black hole data constructed as above, which 
(i) has a conformal factor analogous to Schwarzschild 
at each puncture, and which (ii) has a 1/r fall-off in 
the extrinsic curvature as in (0), we find that the evo- 
lution can also be performed on R 3 , if the metric is 
rescaled with the time-independent conformal factor (f> — 
1 + Yli=i m i/(2ri). We use 1 instead of u so that the 
derivatives of which we need to compute the deriva- 
tives of the unsealed metric g ab , can be computed ana- 
lytically. 

In the numerical implementation, we do not rescale the 
extrinsic curvature but again place the grid such that the 
cl fall between the points of the grid. Note that this is 
unproblematic for the evolution equation (j^) only if the 
shift vector vanishes, because then no derivatives of the 
extrinsic curvature are needed, so let us turn to the choice 
of lapse and shift. 

To complete the specification of the evolution problem, 
we have to choose a slicing condition . Not only can 
lapse and shift be specified freely, but quite non-trivial 
choices have to be made, since simple "Cartesian" choices 
like geodesic slicing (a = 1, (3 a = 0) do not work in 
general. Here we choose vanishing shift, (3 a = 0, and 
maximal slicing for the lapse, where a is the solution of 

Aa = aK ab K ab . (6) 

Below we discuss a few generalizations. 

The question arises whether @ is well-defined near 
the punctures. The principal part of Aa is of order r , 
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while the terms involving first derivatives of a have co- 
efficients of order r 3 , and K ab K ab = i>- 12 8 ac 8 bd K ab K cd 
for puncture data is of order r 6 (r 8 if there is no spin). 
Heuristically, this implies that the first derivatives of a 
have to vanish at the punctures. In the numerics, we 
multiply equation (|^) by <j) 4 ~ 1 /r 4 so that the principal 
part is of order 1, and it turns out that for grids stag- 
gered about the punctures the first derivatives of a are 
sufficiently close to zero near the punctures. 

Factoring out the singular behaviour in the evolution 
equations and in the maximal slicing equation by factor- 
ing out the singular part of the conformal factor of the 
initial data would not be useful if during the course of the 
evolution the character of the singularity at the punctures 
could change or if additional singularities could arise. 
Suppose that a — ► 1 sufficiently fast and that j3 a = 0. 
From (|l|) and (||), near the punctures d t g a b = 0(l/r), and 
also d t K ab = 0(l/r), since -2K ac K c b + K ab K = 0(r 2 ), 
but for g ab = (1 + mj (2r)) i S ab one finds that R rr = 
— &m/(r(m + 2r) 2 ). However, with the time-independent 
rescaling g ab = ip 4 g ab and K ab = i\}^K ab = ip~ 2 K ab , one 
obtains d t g ab = and d t K ab — at the punctures. In 
the numerics, we still rely on the staggering of the grid 
instead of a rescaling of K ab , although the latter may 
have some advantages. 

The bottom line for the puncture method for evolu- 
tions using maximal slicing with vanishing shift is that 
we can treat not only the initial data but also the maxi- 
mal slicing and evolution equations on a domain without 
special inner boundary. As we discuss below, the numer- 
ical implementation shows some remnants due to finite 
difference problems near the punctures, but the general 
setup appears to be valid. Let us add that at the outer 
boundary a Robin boundary condition is used for the ini- 
tial data, for the evolution and maximal slicing the data 
at the outer boundary is kept fixed (cmp. |J). 



III. IMPLEMENTATION AND CODE 
VALIDATION 

The computer code, BAM, is a combination of a 
leap frog evolution code |Tc| ] and a multigrid elliptic solver 
||l8|[ . One of the drawbacks of elliptic slicing conditions 
in 3d codes is that they can account for most of the com- 
putations [0. For our runs, the elliptic solves account for 
60% of the runtime. BAM is outer loop parallelized using 
Power C on an Origin 2000 (distributed shared memory 
model) . A typical run described below can be performed 
on 16 processors in 1-5 GByte of memory in 2-20 hours. 

One of the interesting technical features of BAM is 
that it is the first code used in 3d numerial relativity to 
support fixed mesh refinements. In 3d calculations it is 
an enormous, and perhaps crucial resource savings to, 
for example, use 4 nested cubical boxes centered at the 
origin of width 8, 16, 32, and 64, each with the same 
number of V grid points (4V points total) but with cor- 
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FIG. 1. Convergence results for geodesically sliced 
Schwarzschild at t = 2. The discretized Hamiltonian H(h) 
is shown for both numerical and analytical data. 



respondingly increasing grid spacing, instead of using a 
single big box of size 64 with 128 x AV grid points to ob- 
tain the same resolution at the center. Such fixed mesh 
refinements make it possible to move the outer boundary 
sufficiently far into the 1/r fall-off region. The code is ca- 
pable of (dynamical) adaptive mesh refinement |l(i(| , but 
currently the numerical problems at the grid interfaces 
and problems with parallclization outweigh the potential 
benefits. 

Both the ADM part and the elliptic part of BAM have 
been tested separately for convergence and the propaga- 
tion of the constraints, and compared with analytic or 
linearized results ■ Corresponding tests have been 

performed for the combined code, in particular it repro- 
duces the results on maximal slicing in 3d of M. What 
turned out to be an important consistency and conver- 
gence check, for a given set of parameters the code is 
always run for different resolutions, different refinements 
and different outer boundaries. Since code validation is 
crucial, we will discuss in the remainder of this section 
various issues related to convergence. 

Fig. [l] shows convergence results for the ADM evolution 
in geodesic slicing for a single Schwarzschild black hole 
of mass one. After a finite proper time t = 7r, the initial 
slice reaches the physical singularity. Here we show vari- 
ous plots for the Hamiltonian constraint at t = 2. There 
are two boxes covering the cubical intervals [— 4,4] 3 and 
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[—8, 8] 3 with a refinement factor of 2. Three different cen- 
tral resolutions are considered, 0.250, 0.125, and 0.0625, 
corresponding to 33 3 , 65 3 , and 129 3 points per box, re- 
spectively. 

Second order convergence in the Hamiltonian con- 
straint H is obtained when halving the grid spacing re- 
duces the Hamiltonian constraint by a factor of four, 
H(h) = H(2h)/A. Assuming that a function / can be 
approximated at any grid point and for any grid spacing 
h by f(h) = f + h s e(h), where e(h) is a smooth error 
term that is of order one in h, we can use the standard 
two- and three-level formulas to estimate the non-local 
and local order of convergence, 

s 2 (h)=log 2 \f(2h)\/\f(h)\, (7) 

s 3 (h) = log 2 1/(4/!) - /(2ft) |/|/(2ft) - f(h)\, (8) 

S2(h,x)=log 2 f(2h,x)/f(h,x), (9) 

s 3 (h,x)=log 2 (f(ih,x)-f(2h,x)) 

/(f(2h,x)-f(h,x)), (10) 

where the two-level formulas assume that / is analytically 

zero. We use the ^2-norm \ f{h)\ = (f Eti/C 1 ! 1 ') 2 )^- 

In Fig. 0, we also compare the numerical results with 
the analytically known solution (T^j2^]. Of course, the 
finite difference formula for H(h) computed on analytic 
data is not identically zero, but rather reflects the non- 
vanishing truncation error of the discretization. Note 
that the effect of the fixed outer boundary is clearly vis- 
ible, reducing convergence to zero in a region that actu- 
ally moves inwards at about the speed of light. The inner 
region is causally disconnected from this boundary prob- 
lem and shows rather perfect second order convergence, 
except right at the center, and at points where the log- 
arithm in s 2 becomes ill-defined. Note the small bump 
due to the box nesting interface. 

We now turn to the results obtained with BAM for 
the two black hole problem. Obviously, there is a huge 
parameter space to explore, here we restrict ourselves to 
one particular data set, with the overall scale set by m-i'- 
rm = 1.5, to 2 = 1, 51,2 = (0,0, ±1.5), Pi, 2 = (±2,0,0), 
5i = (-0.5,0.5,0), S 2 = (0,1,1). The initial ADM mass 
estimate is 3.1. 

The basic feature of all maximal slicing runs for black 
hole puncture data is the explosion of the metric and the 
collapse of the lapse, see Fig. || Here we use centered 
cubical grids with 65 3 points each covering [— 4,4] 3 to 
[—32, 32] 3 in spatial coordinates with a constant refine- 
ment factor of 2. The grid spacing therefore ranges from 
0.125 at the center to 1.000 in the outer regions. 

For these parameters, the code runs to about t rj 22. 
For Schwarzschild and widely separated Brill-Lindquist 
holes about t « 25-30 is reached, which is the same 
range as reported for maximally sliced Schwarzschild in 
||. The characteristic feature of maximal slicing is that 
the lapse collapses near the holes, a — > 0, while the slice 
advances with a = 1 at infinity. This has been found 
to avoid the physical singularities for Schwarzschild and 
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FIG. 2. Basic features of the evolution with maximal slic- 
ing and vanishing shift for a general data set: explosion of the 
rescaled metric and collapse of the lapse for t = 0, 2, . . . , 20. 
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FIG. 3. The rescaled metric and the Hamiltonian con- 
straint at t = (general data). 
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FIG. 4. The rescaled metric and the Hamiltonian con- 
straint at t = 10 (general data). 

Misner data, but the "grid-stretching" leads to divergent 
metric coefficients. The code breaks down first when 
solving the maximal slicing equation. The same problem 
occurs at about the same time if a simple SOR (simul- 
taneous overrelaxation) elliptic solver is used instead of 
the multigrid solver, although it is not completely clear 
whether this is a general problem of (|^) or of the numerics 
(cmp. ||). 

Fig. 2 shows that for our data set the lapse collapses 
with roughly the expected speed for the different masses. 
After t « 20, the central region containing both holes 
will not evolve much further, but as we will discuss in 
the next section, the apparent horizon of the two holes 
keeps moving outward. 

In Figs. |3|-|3] we show the rescaled metric and the 
Hamiltonian constraint at times t = 0, 10, and 20. In 
Fig. ||, the three-level convergence based on the norm 
over the z-axis is shown during the course of the run. 
This run uses two grids covering [—8, 8] 3 and [—16, 16] 3 . 
Three different central resolutions are considered, 0.500, 
0.250, and 0.125, corresponding to 33 3 , 65 3 , and 129 3 
points per box, which is chosen a factor of two coarser 
than for the Schwarzschild run in order to move the outer 
boundary outwards by a factor of two. 

At t = 0, Fig. [| the initial data is second order con- 
vergent except right near the punctures. As noted in 
even for rather few points across the punctures, 
the convergence rate away from the punctures is not af- 
fected. When just considering initial data, one can in- 
crease the central resolution and resolve e.g. the rather 
rough humps in the Hamiltonian constraint. For evolu- 
tion problems, both the outer boundary and the inner 
box interfaces become more of a problem, and therefore 
we had to settle on h = 0.125 as the best currently achiev- 
able resolution in the center. 



FIG. 5. The rescaled metric and the Hamiltonian con- 
straint at t = 20 (general data). 

At later times, t = 10 in Fig. | and t = 20 in Fig. 
^, note that the critical region near the punctures does 
not loose convergence, which remains around two for the 
maxima in the Hamiltonian constraint. There is some 
high frequency noise visible in the Hamiltonian constraint 
at t — 20, which is due to the inherent mesh drifting of 
the double leapfrog scheme. For longer run times, it may 
become important to use a more sophisticated scheme. 

The spikes in H at z — ±8 are due to the interpo- 
lations and injections at the box interfaces, which are 
implemented as second order or higher and should lead 
to at least first order in time. Although these spikes do 
no seem to have a profound effect on the evolution, they 
do pose a tricky problem for the following reason, which 
is intrinsic to black hole evolutions with grid-stretching. 
The explosion in the metric always happens faster on 
the coarser, outer grid, so that there is an inherent drift 
between the two grids. This clearly needs further inves- 
tigation. 

Finally, consider Fig. ||. While taking the norm over 
the z-axis smooths out all the local noise, it is still reas- 
suring that there is no catastrophic loss of convergence 
during this run. In particular, the metric converges to at 
least order 1.5 until t = 20. In comparison, for the ana- 
lytic Schwarzschild spacetime discussed earlier, one finds 
convergence of only around 1.8 due to finite size effects. 

Let us mention another problem that effects conver- 
gence. Comparing with geodesic slicing, maximal slicing 
adds the problem that now local errors, in particular at 
the outer boundary and at the punctures, can spread in- 
stantaneously across the grid because of the elliptic char- 
acter of this gauge choice. Even though the trace of K a b 
is roughly second order convergent, Fig. [| we observe 
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FIG. 6. Three-level convergence on the z-axis versus time 
for the rescaled metric, the trace of the extrinsic curvature, 
and the Hamiltonian constraint (general data). 



drifting in the maximal slices, which have effects in the 
time domain similar to the zero crossing artifacts in the 
convergence rate for the Hamiltonian, Fig. []]. 

In summary, based on the above convergence analysis, 
the combined computer code is converging during evolu- 
tions at roughly an order of 1.5 or better. In particular, 
there is no indication that the presence of the punctures 
in the domain of computation destabilizes the code, at 
least for the achievable run times. The run times appear 
to be limited by steep gradients occuring in maximal slic- 
ing away from the punctures. 



IV. RESULTS AND DISCUSSION 

Black hole regions of a spacetime are defined through 
the existence of event horizons, which have been found 
numerically |p3f , but given the available time interval it 
is more promising to study the apparent horizon, which 
is defined for each spatial slice. We therefore look for a 
closed 2d surface S whose outgoing null expansion van- 
ishes, E[S] = 0. Solving such a generalized minimal 
surface equation is an involved problem 24-27]], but for 
our purpose we found that the following simple imple- 
mentation of a curvature flow method p5[ is adequate. 
In terms of the standard Cartesian and spherical coor- 
dinates, a surface which is topologically a sphere can 
be parametrized by u(x,y,z) = r — h(9,ip) = 0, and 
from (20) with surface normal s a — d a u/\du\, \du\ 2 = 
g ab d a udbU, it follows that 



E[u] = (g 



a b 



d a ud b u.. 1 _ 
\ou\ z \ou\ 



K ab ). (11) 



Note that we can evaluate E on a 3d Cartesian grid with- 
out explicitly introducing a 2d surface parametrized by 
(9, ip) since E^ [h] (9, <p) = E^[r - h](r = h,9,ip) = 
E[u](x',y',z') where x' a = ±x a = (1 - f)x a . Starting 
with a large sphere, we iterate u new = u — d\\du\E[u] for 
dX the largest step size that leaves the method stable. 



FIG. 7. Marginally trapped surfaces and g zz in the y — 
plane for t = 20, data from the two innermost nested boxes. 

Fig. [?| shows a 3d visualization of the situation at 
t — 20, while Fig. || sketches the evolution in coordinate 
time t (same data as for Fig. ||). Note that even if the 
maximal t was larger, we cannot expect to see several rev- 
olutions of two "bodies", and that for two reasons: the 
lapse freezes the evolution in the interior, and by con- 
struction the punctures at Ci do not move (see |ll| for 
moving single black holes). 

Since "freezing the evolution of the black holes" is a 
limitation of our method, we want to point out that the 
limitation is far less severe than is immediately apparent. 
The black holes are characterized by a center (the punc- 
ture) and the apparent horizon. The punctures are fixed 
on the grid, but they just correspond to the internal, 
asymptotically flat region of the holes, not the surface of 
the black holes. 

The apparent horizon does move. The lapse collapses 
and freezes the region around the punctures completely, 
but the region near the horizon remains dynamic. In 
head-on collisions of black holes ^ with maximal slicing 
it was found that the gravitational variables are suffi- 
ciently dynamic to allow for the emission of the gravi- 
tational waves that an observer at infinity will measure. 
Given the robustness of wave extraction in 3d demon- 
strated in (ll| for a single distorted black hole and a 
collapsing lapse, one would expect that good wave ex- 
traction is possible for binary black holes in 3d if run 
times of about SOMadm can be achieved, even if one 
uses maximal slicing and vanishing shift. 

Allowing a non-vanishing shift can reduce the strain 
and the shear of the coordinates pl| . A recent pro- 
posal for a minimal strain shift and an algebraic lapse 
to construct approximately corotating coordinates can 
be found in |2q| . Even if the punctures are fixed at a 
coordinate location, one can code the orbital motion in 
the shift vector. Furthermore, it should be possible to 
introduce a translation of the punctures through a shift 
vector. Near the punctures, the conformal factor would 
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FIG. 8. Evolution of the apparent horizon. By construction, the dynamics is contained in the metric while the centers of 
the black holes do not move. Solid lines indicate that a solution to E — has been found, while dotted lines represent surfaces 
which are numerically close to E = everywhere but analytically the flow does not stop and the existence of a true solution 
can be ruled out (based on numerical error estimates for Schwarzschild and Brill-Lindquist data, see also [25]). 
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FIG. 9. Wave indicator log|r 2 p r | along the x-toas, at 
t = 9.25, 12.75, 16.25. A "wave" defined by the zeros in p r 
(— oo in the plot) propagates outwards, passing the numerical 
"noise" near the refinement interface at x — ±8 (cmp. Fig. 5 
of [5]). 



result is shown in Fig. ^ for two nested boxes with 129 3 
points covering [— 8,8] 3 and [— 16, 16] 3 . Note that we do 
not integrate over a sphere to obtain the total flux, which 
would give a much smoother picture. Rather we plot the 
rescaled radial flux on the x-axis to display the local or- 
der of magnitude of the signal and the numerical noise 
at the interior grid faces. In a higher dimensional plot 
of axisymmetric data one finds cubical shaped noise but 
also rather clean axisymmetric signals. Let us emphasize 
that in a binary black hole scenario one can interpret 
the metric in a consistent, gauge invariant manner as a 
wave travelling on a background spacetime only in the 
asymptotically flat region, which is not done here. The 
Bel-Robinson flux is computed only to show that wave- 
like phenomena occur (in particular, they are unrelated 
to the location of the apparent horizon). 

V. CONCLUSION 



not change. Even when the lapse has collapsed to zero, a 
non-vanishing shift leads to motion of the punctures, see 
(0) and (|). _ 

While maximal slicing with its singularity avoing prop- 
erty may be sufficient for some tasks of wave extraction, 
one has to look elsewhere for long time evolutions. A 
simple scenario, although technically involved, is offered 
by apparent horizon boundary conditions. The puncture 
method with maximal slicing can be used to start the 
evolution, while the apparent horizon boundary condi- 
tion takes over at a later time 

To find out whether the code allows wave-like phenom- 
ena, we compute the Bel-Robinson flux, which assuming 
the Einstein equations are satisfied can be expressed as 




p c = E ab e bcd B d a = P ab Q ca \ (12) 

where E ab and B ab axe the electric and magnetic part of 
the conformal tensor, and Q a bc and P a b are the projec- 
tions of the 4d Riemann tensor defined in j2(J . A typical 



In conclusion, we have shown how the standard ap- 
proach of ADM evolution with maximal slicing and van- 
ishing shift can be applied to non-symmetric black hole 
data containing black holes with linear momentum and 
spin by using a time-independent conformal rescaling 
based on the puncture representation of the black holes. 
We discuss an example based on a concrete three dimen- 
sional numerical implementation. The main result of the 
simulations is that this approach allows for the first time 
to evolve through a brief period of the merger phase of the 
black hole inspiral. Looking back to the early numerical 
work of Smarr and Eppley on axisymmetric black hole 
collisions in the seventies ||, we feel that it is useful to 
point out what concretely can be done about the two 
black hole problem today, even if it is just a first step. 

Important issues for further investigations are extend- 
ing the run time so that wave extraction becomes possi- 
ble, non- vanishing shift conditions, and the implementa- 
tion of an apparent horizon boundary condition. These 
issues will be addressed in a new collaborative computa- 
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tional framework, the Cactus code p9|] , which provides 
the infrastructure for input-output and MPI (The Mes- 
sage Passing Interface) based parallelism, and which al- 
lows easy code-sharing among several authors. For ex- 
ample, Cactus includes various evolution modules, initial 
data set constructions, and analysis modules for wave ex- 
traction and apparent horizon finding p7| , p9| . In partic- 
ular, the ADM evolution routines and the multigrid el- 
liptic solver of BAM have already been ported to Cactus, 
adaptive mesh refinement is still under development. An 
important and immediate application of this new, more 
powerful framework will be the comparison of the evolu- 
tion of axisymmetric black hole data in 3d with results 
obtained with 2d codes Q . 

It is a pleasure to thank J. Ehlers, C. Gundlach, P. 
Hiibner, B. Schutz, E. Seidel, P. Walker, and especially S. 
Brandt and B. Schmidt, for their support and many stim- 
ulating discussions. The computations were performed at 
the AEI in Potsdam. 
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